package scu.maqiang.numeric;

import scu.maqiang.csparsej.tdouble.Dcs_common;
import scu.maqiang.csparsej.tdouble.Dcs_entry;
import scu.maqiang.csparsej.tdouble.Dcs_util;
import java.util.*;
import java.util.stream.IntStream;

import static java.lang.Math.*;

import java.io.BufferedWriter;
import java.io.File;
import java.io.FileNotFoundException;
import java.io.FileWriter;
import java.io.IOException;

/**
 * 实数型稀疏矩阵类, 以动态数组的数组形式实现
 * 
 * @author 马强
 * @version 1.0
 * @since 1.0
 *
 */

public class SRMatrix {
	/**
	 * 矩阵行数
	 */
	private int m = 0;
	/**
	 * 矩阵列数
	 */
	private int n = 0;

	/**
	 * 矩阵非0元素个数
	 */
	private int nnz = 0;

	/**
	 * 指标矩阵各行是否行按指标排好序
	 */
	private boolean sorted = false;

	/**
	 * 动态数组的数组
	 */
	public ArrayList<SREntry>[] mat;

	/**
	 * 稀疏矩阵构造函数，设置矩阵行数与列数，矩阵初始为空矩阵
	 * @param m 矩阵行数
	 * @param n 矩阵列数
	 */
	@SuppressWarnings("unchecked")
	public SRMatrix(int m, int n) {
		this.m = m;
		this.n = n;
		this.nnz = 0;
		mat = new ArrayList[m];
		for (int i = 0; i < m; i++) {
			mat[i] = new ArrayList<SREntry>();
		}
	}

	/**
	 * 方形矩阵构造函数，设置方形矩阵阶数
	 * @param n 方形矩阵阶数
	 */
	public SRMatrix(int n) {
		this(n, n);
	}

	/**
	 * 由行列指标构造稀疏矩阵
	 * @param i 矩阵行指标
	 * @param j 矩阵列指标
	 * @param value 矩阵元素值
	 * @param row 矩阵行数
	 * @param col 矩阵列数
	 */
	public SRMatrix(int[] i, int[] j, double[] value, int row, int col) {
		this(row, col);
		for (int ii = 0, nn = i.length; ii < nn; ii++) {
			if (value[ii] != 0.0) {
				mat[i[ii]].add(new SREntry(j[ii], value[ii]));
				this.nnz++;
			}
		}
	}

	/**
	 * 由行列指标构造方形稀疏矩阵
	 * @param i 矩阵行指标
	 * @param j 矩阵列指标
	 * @param value 矩阵元素值
	 * @param size 矩阵阶数
	 */
	public SRMatrix(int[] i, int[] j, double[] value, int size) {
		this(i, j, value, size, size);
	}

	/**
	 * 由已知矩阵构造新矩阵
	 * @param sm 已知矩阵
	 */
	@SuppressWarnings("unchecked")
	public SRMatrix(SRMatrix sm) {
		this(sm.m, sm.n);
		for (int i = 0; i < m; i++) {
			for (SREntry ele : sm.mat[i]) {
				mat[i].add(new SREntry(ele));
				this.nnz++;
			}
		}
	}

	public SRMatrix(double[][] denseMat) {
		this(denseMat.length, denseMat[0].length);
		for(int i = 0; i < m; i++) {
			for(int j = 0; j < n; j++) {
				double ele = denseMat[i][j];
				if (ele != 0.0) {
					mat[i].add(new SREntry(j, ele));
					++nnz;
				}
			}
		}
	}

	/**
	 * 通过外部文件读入稀疏矩阵
	 * 
	 * @param matFile 文件名
	 */
	@SuppressWarnings("unchecked")
	public SRMatrix(String matFile) {
		Scanner sc = null;
		try {
			sc = new Scanner(new File(matFile));
		} catch (FileNotFoundException e) {
			e.printStackTrace();
		}

		int m = sc.nextInt();
		int n = sc.nextInt();
		int nnz = sc.nextInt();
		String isSorted = sc.next();
		sc.next();
		this.m = m;
		this.n = n;
		mat = new ArrayList[m];
		for (int i = 0; i < m; i++) {
			mat[i] = new ArrayList<SREntry>();
		}

		for (int i = 0; i < nnz; i++) {
			int iIdx = sc.nextInt();
			int jIdx = sc.nextInt();
			double value = sc.nextDouble();
			mat[iIdx - 1].add(new SREntry(jIdx - 1, value));
		}

		if (isSorted.trim() == "T") {
			this.sorted = true;
		} else {
			this.sorted = false;
		}
	}

	/**
	 * 拷贝得到一个新矩阵
	 * @return 与原矩阵具有相同结构的新矩阵
	 */
	public SRMatrix copy() {
		return new SRMatrix(this);
	}


	/**
	 * 将稀疏矩阵保存到外部文件中
	 * @param matFile 要保存稀疏矩阵的文件名
	 */
	public void save(String matFile) {
		sort();
		boolean isSymmetric = isSymmetric();
		String isSym = isSymmetric ? "T" : "F";
		try (BufferedWriter bw = new BufferedWriter(new FileWriter(matFile))) {
			bw.write(m + "\t  " + n + "\t  " + nnz + "\t  T\t  " + isSym + "\n");
			for (int i = 0; i < m; i++) {
				for (SREntry ele : mat[i]) {
					bw.write((i + 1) + "\t  " + (ele.col + 1) + "\t  " + ele.val + "\n");
				}
			}
		} catch (IOException e) {
			e.printStackTrace();
		}
	}

	/**
	 * 清除矩阵中所有元素
	 */
	public void clear() {
		for (int i = 0; i < m; i++) {
			mat[i].clear();
			// mat[i] = null;
		}
		nnz = 0;
		sorted = false;
	}
	
	/**
	 * 将矩阵设置为对角矩阵, 并设置对角元素的值
	 * @param val
	 */
	public void setDiag(double val) {
		for (int i = 0; i < m; i++) {
			mat[i].clear();
			mat[i].add(new SREntry(i, val));
		}
		nnz = m;
		sorted = true;
	}

	/**
	 * 交换矩阵两行
	 * @param i 矩阵行指标
	 * @param j 矩阵行指标
	 */
	public void exchangeRow(int i, int j) {
		ParamCheck.checkBounds(i, 0, m);
		ParamCheck.checkBounds(j, 0, n);
		ArrayList<SREntry> temp = mat[i];
		mat[i] = mat[j];
		mat[j] = temp;
	}

	/**
	 * 交换矩阵两列
	 * @param i 矩阵第i列
	 * @param j 矩阵第j列
	 * @param nBand 默认矩阵具有结构, 指示矩阵的半带宽
	 */
	public void exchangeCol(int i, int j, int nBand) {
		ParamCheck.checkBounds(i, 0, m);
		ParamCheck.checkBounds(j, 0, n);
		ArrayList<SREntry> iColTemp = new ArrayList<>();
		ArrayList<SREntry> jColTemp = new ArrayList<>();
		int iIdxBegin = Math.max(0, i - nBand - 1);
		int iIdxEnd = Math.min(m, i + nBand + 1);
		for (int iIdx = iIdxBegin; iIdx < iIdxEnd; iIdx++) {
			for (int jj = 0, n = mat[iIdx].size(); jj < n; jj++) {
				SREntry temp = mat[iIdx].get(jj);
				if (temp.col == j) {
					iColTemp.add(new SREntry(iIdx, temp.val));
					mat[iIdx].remove(jj);
					break;
				}
			}
		}

		// System.out.println(MVO.toString(this.toDenseMatrix()));
		int jIdxBegin = Math.max(0, j - nBand - 1);
		int jIdxEnd = Math.min(m, j + nBand + 1);
		for (int jIdx = jIdxBegin; jIdx < jIdxEnd; jIdx++) {
			for (int ii = 0, n = mat[jIdx].size(); ii < n; ii++) {
				SREntry temp = mat[jIdx].get(ii);
				if (temp.col == i) {
					jColTemp.add(new SREntry(jIdx, temp.val));
					mat[jIdx].remove(ii);
					break;
				}
			}
		}

		// System.out.println(MVO.toString(this.toDenseMatrix()));

		for (SREntry ele : iColTemp) {
			mat[ele.col].add(new SREntry(i, ele.val));
		}

		for (SREntry ele : jColTemp) {
			mat[ele.col].add(new SREntry(j, ele.val));
		}
	}

	/**
	 * 计算第i列元素中从第i个元素以下所有元素的绝对值最大值, 得到该元素值以及对应的行指标
	 * @param i 列数
	 * @param nBand 矩阵带宽
	 * @return 第i列元素中从第i个元素以下所有元素的绝对值最大值, 得到该元素值以及对应的行指标
	 */
	public Pair<Integer, Double> maxAbsColValDown(int i, int nBand) {
		int end = Math.min(i + nBand + 3, m);
		double maxval = getElement(i, i);
		int maxIdx = i;
		for (int ii = i + 1; ii < end; ii++) {
			double temp = getElement(ii, i);
			if (Math.abs(temp) > Math.abs(maxval)) {
				maxval = temp;
				maxIdx = ii;
			}
		}
		return new Pair<Integer, Double>(maxIdx, maxval);
	}

	/**
	 * 将稀疏矩阵转化为一个稠密矩阵
	 * 
	 * @return 二阶数组
	 */
	public double[][] toDenseMatrix() {
		double[][] result = new double[m][n];
		for (int i = 0; i < m; i++) {
			for (SREntry ele : mat[i]) {
				result[i][ele.col] = ele.val;
			}
		}
		return result;
	}

	/**
	 * 得到第i行中列指标为j的SREntry对象
	 * @param i 矩阵第i行
	 * @param j 矩阵第j列
	 * @return 返回SREntry对象, 对象中col分量为j, val为元素的值, 
	 *         如果稀疏矩阵中不存在第i行第j列元素, 则返回null
	 */
	public SREntry getSREntry(int i, int j) {
		ParamCheck.checkBounds(i, 0, m);
		ParamCheck.checkBounds(j, 0, n);
		for (SREntry ele : mat[i]) {
			if (ele.col == j) {
				return ele;
			}
		}
		return null;
	}

	/**
	 * 得到矩阵元素的值
	 * @param i 元素行指标
	 * @param j 元素列指标
	 * @return 如果存在，则返回元素的值；否则返回0
	 */
	public double getElement(int i, int j) {
		ParamCheck.checkBounds(i, 0, m);
		ParamCheck.checkBounds(j, 0, n);
		Optional<SREntry> element = mat[i].stream().filter(ele -> ele.col == j).findFirst();
		return element.isPresent() ? element.get().val : 0.0;
	}

	/**
	 * 设置矩阵元素的值
	 * @param i 元素行指标
	 * @param j 元素列指标
	 * @param value 元素值
	 * @return 当前稀疏矩阵
	 */
	public SRMatrix setElement(int i, int j, double value) {
		ParamCheck.checkBounds(i, 0, m);
		ParamCheck.checkBounds(j, 0, n);
		int size = mat[i].size();
		if (size == 0 && value != 0.0) {
			mat[i].add(new SREntry(j, value));
			this.nnz++;
			return this;
		}

		for (int idx = 0; idx < size; idx++) {
			SREntry se = mat[i].get(idx);
			if (se.col == j) {
				if (value != 0.0)
					se.val = value;
				else {
					mat[i].remove(idx);
					--(this.nnz);
				}
				return this;
			}
		}

		mat[i].add(new SREntry(j, value));
		++(this.nnz);
		return this;
	}

	/**
	 * 不作任何检查设置矩阵元素的值
	 * @param i 元素行指标
	 * @param j 元素列指标
	 * @param value 元素值
	 * @return 当前稀疏矩阵
	 */
	public SRMatrix setElementFast(int i, int j, double value) {
		ParamCheck.checkBounds(i, 0, m);
		ParamCheck.checkBounds(j, 0, n);
		mat[i].add(new SREntry(j, value));
		++(this.nnz);
		return this;
	}

	/**
	 * 在稀疏矩阵相应行列指标上加上相应元素值, 最终结果存在三种情况:
	 * 1. 如果元素值为0，则不做任何操作
	 * 2. 如果元素值非0，并且矩阵中行列指标中存在元素，则将原来元素值与之相加,
	 *    所加和如果为0, 则删除该元素，如果非0，则修改该元素值
	 * 3. 如果元素值非0，但矩阵中不存在该元素，则向矩阵中新增加该元素值
	 * @param i 元素行指标
	 * @param j 元素列指标
	 * @param val 元素值
	 *            
	 */
	public void addElement(int i, int j, double val) {
		ParamCheck.checkBounds(i, 0, m);
		ParamCheck.checkBounds(j, 0, n);
		if (val != 0.0) {
			int size = mat[i].size();
			if (size == 0) {
				mat[i].add(new SREntry(j, val));
				this.nnz++;
				return;
			}
			for (int idx = 0; idx < size; idx++) {
				SREntry se = mat[i].get(idx);
				if (se.col == j) {
					double temp = se.val + val;
					double minval = Math.min(Math.abs(se.val), Math.abs(val));
					double ratio = Math.abs(temp) / minval;
					if (ratio <  2.47009e-10) {
						//System.out.println(se.val + "  " + val);
						//System.out.println("OK!");
						mat[i].remove(idx);
						--(this.nnz);
					} else {
						se.val = temp;
					}
					return;
				}
			}
			mat[i].add(new SREntry(j, val));
			++(this.nnz);
		}
	}

	/**
	 * 在稀疏矩阵相应行列指标上加上相应元素值, 不检查加上后的结果为0的情况
	 * @param i
	 * @param j
	 * @param val
	 */
	public void addElementFast(int i, int j, double val) {
		ParamCheck.checkBounds(i, 0, m);
		ParamCheck.checkBounds(j, 0, n);
		int size = mat[i].size();
		if (size == 0) {
			mat[i].add(new SREntry(j, val));
			this.nnz++;
			return;
		}
		for (int idx = 0; idx < size; idx++) {
			SREntry se = mat[i].get(idx);
			if (se.col == j) {
				se.val += val;
				return;
			}
		}
		mat[i].add(new SREntry(j, val));
		++(this.nnz);
	}

	/**
	 * 在矩阵中删除第i行第j列所有元素值，这里假设稀疏矩阵中所有非0元素的行列指标是对称的
	 * @param i 行指标
	 */
	public void removeSymmetricIJ(int i) {
		mat[i].forEach(ele -> {
			int col = ele.col;
			if (col != i) {
				removeElement(col, i);
			}
		});
		removeRow(i);
	}
	
	public void removeRowCol(int i) {
		SRMatrix AT = this.transpose();
		removeRow(i);
		AT.mat[i].forEach(ele -> {
			int row = ele.col;
			removeElement(row, i);
		});
	}

	/**
	 * 去除稀疏矩阵中某行列指标下的元素
	 * @param i 元素行指标
	 * @param j 元素列指标
	 * @return 如果删除的元素存在，则返回true; 否则返回false
	 */
	public boolean removeElement(int i, int j) {
		ParamCheck.checkBounds(i, 0, m);
		ParamCheck.checkBounds(j, 0, n);
		int size = mat[i].size();
		for (int idx = 0; idx < size; idx++) {
			if (mat[i].get(idx).col == j) {
				mat[i].remove(idx);
				--(this.nnz);
				return true;
			}
		}
		return false;
	}

	/**
	 * 删除稀疏矩阵某行所有元素
	 * @param i 矩阵行指标
	 */
	public void removeRow(int i) {
		ParamCheck.checkBounds(i, 0, m);
		this.nnz -= mat[i].size();
		mat[i].clear();
	}

	/**
	 * 得到矩阵行数
	 * @return 矩阵行数
	 */
	public int getM() {
		return m;
	}

	/**
	 * 得到矩阵列数
	 * 
	 * @return 矩阵列数
	 */
	public int getN() {
		return n;
	}

	/**
	 * 得到稀疏矩阵中非0元素个数
	 * 
	 * @return 非0元素个数
	 */
	public int getNNZ() {
		return nnz;
	}

	/**
	 * 得到稀疏矩阵中下三角元素中非0元素个数
	 * @return 下三角非0元素个数
	 */
	public int getNNZL() {
		int sum = 0;
		for (int i = 0; i < m; i++) {
			int ii = i;
			sum += mat[i].stream().filter(ele -> ele.col < ii).count();
		}
		return sum;
	}

	/**
	 * 得到稀疏矩阵中上三角元素中非0元素个数
	 * @return 上三角非0元素个数
	 */
	public int getNNZU() {
		int sum = 0;
		for (int i = 0; i < m; i++) {
			int ii = i;
			sum += mat[i].stream().filter(ele -> ele.col > ii).count();
		}
		return sum;
	}

	/**
	 * 分别得到稀疏矩阵中上三角与下三角非0元素的个数
	 * @return 整数对, 第一个值为上三角非0元素个数, 第二个值为下三角非0元素个数
	 */
	public Pair<Integer, Integer> getNNZLU() {
		int numL = 0;
		int numU = 0;
		for (int i = 0; i < m; i++) {
			int ii = i;
			numL += mat[i].stream().filter(ele -> ele.col < ii).count();
			numU += mat[i].stream().filter(ele -> ele.col < ii).count();
		}
		return new Pair<Integer, Integer>(numL, numU);
	}

	/**
	 * 得到矩阵第i行的动态数组
	 * @param i 矩阵行指标
	 * @return 第i行的动态数组
	 */
	public ArrayList<SREntry> getRow(int i) {
		ParamCheck.checkBounds(i, 0, m);
		return mat[i];
	}

	/**
	 * 实现稀疏矩阵相加, 如果当前矩阵为A, 输入参数为alpha与B, 最终得返回矩阵值C为
	 *    C ＝ A + alpha * B.
	 * @param alpha  实数
	 * @param B　相加矩阵
	 * @return　返回矩阵值为C = A + alpha * B
	 */
	public SRMatrix add(double alpha, SRMatrix B) {
		SRMatrix result = new SRMatrix(this);
		int M = result.m;
		boolean flag = false;
		if (alpha == 0.0) {
			return result;
		}

		for (int i = 0; i < M; i++) {
			int mat2Size = B.mat[i].size();
			for (int j = 0; j < mat2Size; j++) {
				// flag记录是否在第二个加数中存在该指标
				flag = false;
				SREntry e2 = B.mat[i].get(j);
				int resultSize = result.mat[i].size();
				for (int k = 0; k < resultSize; k++) {
					SREntry e1 = result.mat[i].get(k);
					// 寻找两矩阵中的相同行中是否有相同的列指标
					if (e1.col == e2.col) {
						double temp2 = e1.val + alpha * e2.val;
						// 如果有相同列指标，两者之和相加又不为0,则修改列指标下值ֵ
						if (temp2 != 0.0) {
							e1.val = temp2;
						} else {
							//如果相加为0,则删除该元素
							result.mat[i].remove(k);
							result.nnz--;
						}
						flag = true;
						break;
					}
				}
				//如果没有该指标，则直接加入该元素
				if (!flag) {
					result.mat[i].add(new SREntry(e2.col, alpha * e2.val));
					result.nnz++;
				}
			}
		}
		return result;
	}

	/**
	 * 实现稀疏矩阵相加, 如果当前矩阵为A, 输入参数为alpha与B, 最终得返回矩阵值C为
	 *    C ＝ alpha * A + alpha * B.
	 * @param alpha  实数 要求不为0
	 * @param A 相加矩阵
	 * @param beta 实数 要求不为0
	 * @param B　相加矩阵
	 * @return　返回矩阵值为C = alpha * A + beta * B
	 */
	public static SRMatrix Add(double alpha, SRMatrix A, double beta, SRMatrix B) {
		assert(alpha != 0.0);
		assert(beta != 0.0);
		SRMatrix result = new SRMatrix(A.m, A.n);
		int m = A.m;
		for (int i = 0; i < m; i++) {
			for (SREntry ele : A.mat[i]) {
				result.mat[i].add(new SREntry(ele.col, alpha * ele.val));
				result.nnz++;
			}
		}
		boolean flag = false;
		for (int i = 0; i < m; i++) {
			int mat2Size = B.mat[i].size();
			for (int j = 0; j < mat2Size; j++) {
				flag = false;
				SREntry e2 = B.mat[i].get(j);
				int resultSize = result.mat[i].size();
				for (int k = 0; k < resultSize; k++) {
					SREntry e1 = result.mat[i].get(k);
					if (e1.col == e2.col) {
						double temp2 = e1.val + beta * e2.val;
						if (temp2 != 0.0) {
							e1.val = temp2;
						} else {
							result.mat[i].remove(k);
							result.nnz--;
						}
						flag = true;
						break;
					}
				}
				if (!flag) {
					result.mat[i].add(new SREntry(e2.col, beta * e2.val));
					result.nnz++;
				}
			}
		}
		return result;		
	}

	public int[][] getRowRange() {
		int[][] rowRange = new int[m][2];
		sort();
		for(int i = 0; i < m; i++) {
			int size = mat[i].size();
			if (size > 0) {
				rowRange[i][0] = mat[i].get(0).col;
				rowRange[i][1] = mat[i].get(size - 1).col;
			} else {
				rowRange[i][0] = -1;
			}

		}
		return rowRange;
	}

	public int[][] getColRange() {
		int[][] colRange = new int[n][];
		for(int i = 0; i < n; i++) {
			colRange[i] = new int[]{m + 1, -1};
		}
		sort();
		for(int i = 0; i < m; i++) {
			for(SREntry ele: mat[i]) {
				if(colRange[ele.col][0] > i) {
					colRange[ele.col][0] = i;
				}
				if(colRange[ele.col][1] < i) {
					colRange[ele.col][1] = i;
				}
			}
		}
		return colRange;
	}

	public SRMatrix mul(SRMatrix b) {
		assert(n == b.m);
		int bn = b.n;
		SRMatrix result = new SRMatrix(m, b.n);
		for(int i = 0; i < m; i++) {
			for(SREntry eleA: mat[i]) {
				for(SREntry eleB: b.mat[eleA.col]) {
					result.addElement(i, eleB.col, eleA.val * eleB.val);
				}
			}
		}
		return result;
	}

	public SCMatrix mul(SCMatrix b) {
		assert(n == b.getM());
		int bn = b.getN();
		SCMatrix result = new SCMatrix(m, bn);
		Complex[] temp = Complex.zeros(bn);
		for(int i = 0; i < m; i++) {
			Complex.fill(temp, 0.0);
			for(SREntry eleA: mat[i]) {
				for(SCEntry eleB: b.mat[eleA.col]) {
					temp[eleB.col].real += eleA.val * eleB.val.real;
					temp[eleB.col].image += eleA.val * eleB.val.image;
				}
			}
			for(int k = 0; k < bn; k++) {
				if (temp[k].norm() != 0.0) {
					result.setElementFast(i, k, new Complex(temp[k].real, temp[k].image));
				}
			}
		}
		return result;
	}

	public void mul(SRMatrix b, SRMatrix result) {
		assert(m == result.m);
		assert(n == b.m);
		assert(b.n == result.n);
		int bn = b.n;
		int[][] rowRange = getRowRange();
		int[][] colRange = b.getColRange();
		double sum = 0;
		for(int i = 0; i < m; i++) {
			int rowMin = rowRange[i][0];
			int rowMax = rowRange[i][1];
			for(int j = 0; j < bn; j++) {
				sum = 0.0;
				int colMin = colRange[j][0];
				int colMax = colRange[j][1];
				if (rowMax < colMin || colMax < rowMin) {
					continue;
				} else {
					for(SREntry ele: mat[i]) {
						sum += ele.val * b.getElement(ele.col, j);
					}
					if (sum != 0.0) {
						result.setElementFast(i, j, sum);
					}
				}
			}
		}
	}

	/**
	 *  矩阵数量乘法
	 * @param d 实数
	 * @return 返回乘以d后得到的矩阵
	 */
	public SRMatrix mul(double d) {
		assert(d!= 0.0);
		SRMatrix result = new SRMatrix(m, n);
		mmul(d, result);
		return result;
	}

	/**
	 *
	 * @param d
	 * @param result
	 */
	public void mul(double d, SRMatrix result) {
		ParamCheck.checkEqual(result.m, m);
		ParamCheck.checkEqual(result.n, n);
		mmul(d, result);
	}

	private void mmul(double d, SRMatrix result) {
		IntStream.range(0, m).forEach(i -> mat[i].forEach(ele -> {
			result.mat[i].add(new SREntry(ele.col, ele.val * d));
			result.nnz++;
		}));
		/*
		 * for (int i = 0; i < m; i++) { for (SREntry ele : mat[i]) {
		 * result.mat[i].add(new SREntry(ele.col, ele.val * d)); result.nnz++; } }
		 */
	}

	public double[] mul(double[] u) {
		ParamCheck.checkEqual(n, u.length);
		double[] v = new double[m];
		mmul(u, v);
		return v;
	}

	public void mul(double[] u, double[] v) {
		ParamCheck.checkEqual(n, u.length);
		ParamCheck.checkEqual(m, v.length);
		mmul(u, v);
	}

	private void mmul(double[] u, double[] v) {
		for (int i = 0; i < m; i++) {
			v[i] = 0.0;
			for (SREntry ele : mat[i]) {
				v[i] += ele.val * u[ele.col];
			}
		}
	}

	public double[] mul(double coef, double[] u) {
		ParamCheck.checkEqual(n, u.length);
		double[] v = new double[m];
		mmul(coef, u, v);
		return v;
	}

	public void mul(double coef, double[] u, double[] v) {
		ParamCheck.checkEqual(n, u.length);
		ParamCheck.checkEqual(m, v.length);
		mmul(coef, u, v);
	}

	private void mmul(double coef, double[] u, double[] v) {
		// Arrays.parallelSetAll(v, i -> mat[i].parallelStream().mapToDouble(ele ->
		// ele.val * u[ele.col]).sum());
		for (int i = 0; i < m; i++) {
			v[i] = 0.0;
			for (SREntry ele : mat[i]) {
				v[i] += ele.val * u[ele.col] * coef;
			}
		}
	}

	public void mul(Complex[] u, Complex[] v) {
		ParamCheck.checkEqual(n, u.length);
		ParamCheck.checkEqual(m, v.length);
		mmul(u, v);
	}

	private void mmul(Complex[] u, Complex[] v) {
		// Arrays.parallelSetAll(v, i -> mat[i].parallelStream().mapToDouble(ele ->
		// ele.val * u[ele.col]).sum());
		for (int i = 0; i < m; i++) {
			v[i].real = 0.0;
			v[i].image = 0.0;
			for (SREntry ele : mat[i]) {
				v[i].real += ele.val * u[ele.col].real;
				v[i].image += ele.val * u[ele.col].image;
			}
		}
	}

	public double[] mulT(double[] u) {
		ParamCheck.checkEqual(m, u.length);
		double[] v = new double[n];
		mmulT(u, v);
		return v;
	}

	public void mulT(double[] u, double[] v) {
		ParamCheck.checkEqual(m, u.length);
		ParamCheck.checkEqual(n, v.length);
		mmulT(u, v);
	}

	/**
	 * 
	 * @param u
	 * @param v
	 */
	private void mmulT(double[] u, double[] v) {
		Arrays.fill(v, 0.0);
		for (int i = 0; i < m; i++) {
			for (SREntry ele : mat[i]) {
				v[ele.col] += ele.val * u[i];
			}
		}
	}
	
	/**
	 * 
	 * @param rowDof
	 * @param colDof
	 * @param ele
	 */
	public void assemble(int[] rowDof, int[] colDof, double[][] ele) {
		int m = ele.length;
		int n = ele[0].length;
		ParamCheck.checkEqual(m, rowDof.length);
		ParamCheck.checkEqual(n, colDof.length);
		for (int i = 0; i < m; i++) {
			for (int j = 0; j < n; j++) {
				addElement(rowDof[i], colDof[j], ele[i][j]);
			}
		}

	}
	
	/**
	 * 
	 * @param dof
	 * @param ele
	 */
	public void assemble(int[] dof, double[][] ele) {
		assemble(dof, dof, ele);
	}

	public void assemble(int row, int[] colDof, double[] ele) {
		int m = ele.length;
		ParamCheck.checkEqual(m, colDof.length);
		for (int i = 0; i < m; i++) {
			addElement(row, colDof[i], ele[i]);
		}
	}

	/**
	 * 
	 * @param rowDof
	 * @param col
	 * @param ele
	 */
	public void assemble(int[] rowDof, int col, double[] ele) {
		int m = ele.length;
		ParamCheck.checkEqual(m, rowDof.length);
		for (int i = 0; i < m; i++) {
			addElement(rowDof[i], col, ele[i]);
		}
	}

	/**
	 * 
	 */
	public void sort() {
		// for (int i = 0; i < m; i++) {
		// Collections.sort(mat[i], (a, b) -> a.col > b.col ? 1 : -1);
		// }
		IntStream.range(0, m).parallel().forEach(i -> Collections.sort(mat[i], (a, b) -> a.col > b.col ? 1 : -1));
		sorted = true;
	}

	/**
	 * 
	 * @return
	 */
	public boolean isSorted() {
		return sorted;
	}

	/**
	 * 
	 * @return
	 */
	public boolean makeEntrySymmetric() {
		int oldNNZ = nnz;
		for (int i = 0; i < m; i++) {
			int len = mat[i].size();
			for (int j = 0; j < len; j++) {
				int col = mat[i].get(j).col;
				if (getElement(col, i) == 0.0) {
					setElementFast(col, i, 0.0);
				}
			}
		}
		return oldNNZ != nnz;
	}

	@Override
	public String toString() {
		StringBuilder sb = new StringBuilder();
		sb.append(m).append("\t  ").append(n).append("\t  ").append(nnz).append("\t  ").append(isSymmetric()).append("\n");
		for (int i = 0; i < m; i++) {
			for (SREntry ele : mat[i]) {
				sb.append(i + 1).append("\t  ").append(ele.col + 1).append("\t  ").append(ele.val).append("\n");
			}

		}
		return sb.toString();
	}

	public String displayByRow() {
		StringBuilder sb = new StringBuilder();
		sb.append(m + "\t  " + n + "\t  " + nnz + "\t  " + isSymmetric() + "\n");
		for (int i = 0; i < m; i++) {
			for (SREntry ele : mat[i]) {
				sb.append(i + "\t  " + ele.col + "\t  " + ele.val + "\n");
			}

		}
		return sb.toString();
	}

	public String displayByCol() {
		StringBuilder sb = new StringBuilder();
		sb.append(m + "\t" + n + "\t" + nnz + "\t" + isSymmetric() + "\n");
		ArrayList<Pair<Integer, Double>>[] cscmat = new ArrayList[n];
		for(int i = 0; i < n; i++) {
			cscmat[i] = new ArrayList<>();
		}

		for(int i = 0; i < m; i++) {
			for(SREntry ele: mat[i]) {
				cscmat[ele.col].add(new Pair<>(i, ele.val));
			}
		}
		IntStream.range(0, n).parallel().forEach(i -> Collections.sort(cscmat[i], (a, b) -> a.first > b.first ? 1 : -1));
		for(int i = 0; i < n; i++) {
			for(Pair<Integer, Double> ele: cscmat[i]) {
				sb.append(ele.first + "\t  " + i + "\t  " + ele.second + "\n");
			}
		}
		return sb.toString();
	}

	public double[] diag() {
		double[] d = new double[m];
		for (int i = 0; i < m; i++) {
			for (SREntry ele : mat[i]) {
				if (ele.col == i) {
					d[i] = ele.val;
				}
			}
		}
		return d;
	}

	public SRMatrix transpose() {
		SRMatrix result = new SRMatrix(n, m);
		for (int i = 0; i < m; i++) {
			for (SREntry ele : mat[i]) {
				result.mat[ele.col].add(new SREntry(i, ele.val));
				result.nnz++;
			}
		}
		return result;
	}

	public SRMatrix triU() {
		SRMatrix upper = new SRMatrix(m, n);
		for (int i = 0; i < m; i++) {
			int ii = i;
			mat[i].stream().filter(ele -> ele.col > ii).forEach(ele -> upper.setElementFast(ii, ele.col, ele.val));
		}
		return upper;
	}

	public SRMatrix triL() {
		SRMatrix lower = new SRMatrix(m, n);
		for (int i = 0; i < m; i++) {
			int ii = i;
			mat[i].stream().filter(ele -> ele.col < ii).forEach(ele -> lower.setElementFast(ii, ele.col, ele.val));
			// lower.mat[i].addAll(mat[i].stream().filter(ele -> ele.col < ii));
		}
		return lower;
	}

	/**
	 * 使用三数组得到矩阵各元素的值
	 * @param ii　矩阵行指标
	 * @param jj　矩阵列指标
	 * @param values　矩阵元素值
	 */
	public void getTriValues(int[] ii, int[] jj, double[] values) {
		int in = 0;
		for (int i = 0; i < m; i++) {
			for (SREntry ele : mat[i]) {
				ii[in] = i;
				jj[in] = ele.col;
				values[in] = ele.val;
				in++;
			}
		}
	}

	/**
	 * 计算得到稀疏矩阵元素最大值
	 * @return　稀疏矩阵元素最大值
	 */
	public double maxElement() {
		double[] maxval = new double[m];
		for (int i = 0; i < m; i++) {
			maxval[i] = mat[i].stream().mapToDouble(ele -> ele.val).max().orElse(0.0);
		}
		return Arrays.stream(maxval).max().orElse(0.0);
	}

	/**
	 * 计算得到稀疏矩阵元素最大值，并返回最大值所在行列值
	 * @param IJ　最大值所在行列值
	 * @return　稀疏矩阵元素最大值
	 */
	public double maxElement(int[] IJ) {
		double[] maxval = new double[m];
		int[] maxJ = new int[m];
		for (int i = 0; i < m; i++) {
			if (mat[i] != null) {
				int tempJ = mat[i].get(0).col;
				double tempVal = mat[i].get(0).val;
				for (int j = 1, size = mat[i].size(); j < size; j++) {
					SREntry ele = mat[i].get(j);
					if (ele.val > tempVal) {
						tempJ = ele.col;
						tempVal = ele.val;
					}
				}
				maxJ[i] = tempJ;
				maxval[i] = tempVal;
			}
		}
		double result = maxval[0];
		IJ[0] = 0;
		IJ[1] = maxJ[0];
		for (int i = 1; i < m; i++) {
			if (result < maxval[i]) {
				result = maxval[i];
				IJ[0] = i;
				IJ[1] = maxJ[i];
			}
		}
		return result;
	}

	public boolean isSymmetric() {
		for (int i = 0; i < m; i++) {
			for (SREntry ele : mat[i]) {
				if (abs(ele.val - getElement(ele.col, i)) > Constants.Er) {
					return false;
				}
			}
		}
		return true;
	}

	public int[] symamd() {
//			if(isSymmetric()) {
//				int[] degree = new int[m];
//				Arrays.setAll(degree, i -> mat[i].size());
//				System.out.println(Arrays.toString(degree));
//				int[] perm = new int[m];
//				Algorithms.selectionSort(degree, perm, true);
//				System.out.println(Arrays.toString(perm));
//			} else {
//				System.out.println("The matrix is not symmetric");
//	//			return null;
//			}
//		int[] degree = new int[m];
//		boolean[] flag = new boolean[m];
//		Arrays.setAll(degree, i -> mat[i].size());
//		Pair<Integer, Integer> info = new Pair<>();
//		ArrayList<Integer> idxList = new ArrayList<>();
//		for(int i = 0; i < m; i++) {
//			if (flag[i] == false) {
//				for(SREntry ele : mat[i]) {
//					if (flag[ele.col] == false) {
//						
//					}
//				}
//			}
//			//for(int j = i; j < m; j++) {
//			//	
//			//}
//		}
//		
//		Arrays.setAll(degree, i -> mat[i].size());
//		Pair<Integer, Integer> info = new Pair<>();
//		
//		MVO.min(degree, info);
//		int idx = info.getFirst();
//		for(SREntry ele: mat[idx]) {
//			flag[ele.col] = true;
//		}

		int[] degree = new int[m];
		int[] perm = new int[m];
		boolean[] flag = new boolean[m];
		Arrays.setAll(degree, i -> mat[i].size());
		Pair<Integer, Integer> info = new Pair<>();
		for (int i = 0; i < m; i++) {
			MVO.min(degree, info);
			int idx = perm[i] = info.getFirst();
			for (SREntry ele : mat[idx]) {
				if (flag[ele.col] == false) {
					degree[ele.col]--;
				}
			}
			flag[idx] = true;
			degree[idx] = Integer.MAX_VALUE;
		}
		/*
		 * System.out.println(Arrays.toString(perm)); flag = new boolean[m]; for(int i =
		 * 0; i < m; i++) { if (flag[i] == false && i != perm[i]) {
		 * System.out.println("change: " + i + "\t, " + perm[i]); this.exchangeRow(i,
		 * perm[i]); System.out.println(MVO.toString(this.toDenseMatrix()));
		 * this.exchangeCol(i, perm[i], nBand);
		 * System.out.println(MVO.toString(this.toDenseMatrix()));
		 * System.out.println("=================="); flag[i] = true; flag[perm[i]] =
		 * true; } }
		 */

		// int[] inversePerm = new int[m];
		// Arrays.setAll(inversePerm, i -> perm[i]);

		// System.out.println(Arrays.toString(perm));
		SRMatrix temp = new SRMatrix(m);
		for (int i = 0; i < m; i++) {
			for (SREntry ele : mat[perm[i]]) {
				temp.mat[i].add(new SREntry(ele.col, ele.val));
			}
		}

		// System.out.println(MVO.toString(temp.toDenseMatrix()));
		int[] inversePerm = new int[m];
		for (int i = 0; i < m; i++) {
			inversePerm[perm[i]] = i;
		}
		// System.out.println(Arrays.toString(inversePerm));
		this.clear();
		for (int i = 0; i < m; i++) {
			for (SREntry ele : temp.mat[i]) {
				mat[i].add(new SREntry(inversePerm[ele.col], ele.val));
			}
		}
		temp.clear();
		return perm;
	}

	public void scaleRow(double[] ratio) {
		assert(m == ratio.length);
		for(int i = 0; i < m; i++) {
			for(SREntry ele: mat[i]) {
				ele.val *= ratio[i];
			}
		}
	}

	public void scaleCol(double[] ratio) {
		assert(n == ratio.length);
		for(int i = 0; i < m; i++) {
			for(SREntry ele: mat[i]) {
				ele.val *= ratio[ele.col];
			}
		}
	}

	public Dcs_common.Dcs toDcs_matrix() {
		Dcs_common.Dcs T;
		T = Dcs_util.cs_spalloc(0, 0, 1, true, true); /* allocate result */
		String line;
		for(int k = 0; k < m; k++) {
			for(SREntry ele: mat[k]) {
				if (!Dcs_entry.cs_entry(T, k, ele.col, ele.val))
					return (null);
			}

		}
		return (T);
	}

	/**
	 * 实际稀疏矩阵的不完全cholesky分解
	 * @return
	 */
	public SRMatrix iCholesky() {
		assert(m == n);
		SRMatrix R = new SRMatrix(n);

		for(int i = 0; i < m; i++) {
			double temp = 0.0;
			for(int j = 0; j < i; j++) {
				double Rji = R.getElement(j, i);
				temp += Rji * Rji;
			}
			double temp2 = getElement(i, i) - temp;
			if (temp2 <= 0) {
				throw new IllegalArgumentException("THe matrix is not positive definite");
			}
			temp2 = Math.sqrt(temp2);
			R.setElementFast(i, i, temp2);
			for(SREntry ele : mat[i]) {
				if (ele.col > i) {
					temp = 0.0;
					for(int k = 0; k < i; k++) {
						temp += getElement(k, i) * getElement(k, ele.col);
					}
					R.setElementFast(i, ele.col, (ele.val - temp) / temp2);
				}
			}
		}
		return R;
	}
}
